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A method for the construction and implementation of a mathematical zonal model of heat exchange in a 
glass-melting furnace with a horseshoe-shaped flame and an output of 300 ton/day is considered. The boun¬ 
dary conditions of modeling are described in detail. The authors give the results of the computation of external 
heat exchange and hydrodynamics of the melting tank, as well as regression equations of temperature fields 
for refractory brickwork, the flame space, and the glass melt surface. 


The practical implementation of the mathematical model 
of a glass-melting furnace, whose main principles are de¬ 
scribed in [1], implies the testing of the calculation proce¬ 
dures and the fonnation of a database of boundary conditions 
corresponding to the actual parameters of industrial analogs. 
The furnace selected for testing the methodological algo¬ 
rithm and adapting the mathematical model is the regenera¬ 
tor glass-melting furnace with a horseshoe-shaped flame 
with an output of 300 ton/day designed by V. Ya. Dzyuzer. 
The furnace is intended for melting container glass (KT-1, 
ZT-1, BT-1, etc.) with the option of modifying glass melt tint 
in the course of its operation. Below are certain specifica¬ 
tions of this furnace used as boundary conditions in the de¬ 
velopment of its mathematical model. 

The surface area of the melting tank is 115.77 m 2 (length 
and width are 13.62 and 8.50 m, respectively). Fuel: natural 
gas with lower heat-generating capacity = 35.042 MJ/m 3 . 

The rated thermal load of the furnace for glass ZT-1 for air 
heating temperature 1200°C and the batch : cullet ratio equal 
to 70 : 30 is BQ l p = 16.867 MW. Thus, for fuel consumption 

B = 1732.9 m 3 /h the specific heat rate and specific glass melt 
output are 1159.4 kcal/kg and 2.591 ton/(m 2 • day), respec¬ 
tively. 

The concept of the mathematical model of the glass¬ 
melting furnace [1] implies a consecutive solution of the ex¬ 
terior and interior problems. Solving the problem of exterior 
heat exchange in the working space includes several stages. 

The first stage is the construction of a zonal model of the 
furnace and the calculation of the generalized angular radia¬ 
tion coefficient \|/.. (/ is the number of the radiating zone and 
i is the number of the recipient zone). The working space of 
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the furnace and the enclosing brickwork and tank surfaces 
are split into volumetric (gaseous) and surface zones. The 
quantity and sizes of zones are selected to register as pre¬ 
cisely as possible variations in heat and mass exchange pa¬ 
rameters, the tank hydrodynamics, and the stages of the tech¬ 
nological process distributed along the melting tank. 

It our case the working space of the furnace is repre¬ 
sented in the form of a parallelepiped, whose length and 
width correspond to the real sizes of the industrial analog. 
The height is selected based on the condition of the equality 
between the surface areas of the brickwork and the tank: 
F b /F t = const for the furnace and for the model. The work¬ 
ing space of the furnace is spilt into five zones along the 
tank, five zones across its width, and three zones across its 
height (Fig. la). The zone coordinates are as follows: along 
the furnace length X Q ->X 5 = 0, 2.043, 4.086, 7.264, 10.442, 
and 13.620 m; in width Y 0 -A Y 5 = 0, 1.15, 3.35, 5.15, 7.35, 
and 8.50 m, and in height Z Q —> Z 5 = 0, 0.85, 1.70, and 
2.55 m. The limiting surface of the model from below is the 
glass melt surface. 

Thus, the working space of the furnace, including the 
glass melt surface is represented as 75 volumetric and 110 
surface zones. The latter are subdivided into surface zones of 
types I, II, and III. The surface zones of type I are surfaces 
with a preset temperatures. These are the surfaces through 
which air (1350°C) and batch (1450°C) are fed into the fur¬ 
nace and through which the combustion products (1500°C) 
are removed. The glass melt surface represents surface zones 
of type II, for which a distribution of heat flows is specified. 
All enclosing brickwork surfaces are zones of type Ill, for 
which a relation between heat flows and temperature is de¬ 
termined. 
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Fig. 1 . Approximated zonal model of the working space of the 
glass-melting furnace with horseshoe flame direction (a) and longi¬ 
tudinal section of the melting tank ( b). 


The length /, and the width of the melting tank are given 
above. For tinted glass the depth of the tank before and after 
the overflow edge is h { = 1.3 and i 2 =1.6m, respectively 
(Fig. 1 b). The overflow edge of height /t 3 = 0.8 and width 
/ 3 = 0.4 m is placed at a distance of l 2 = 9.2 m from the 
burner inlet. The glass melt moves from the level of the tank 
bottom through a neck of width 0.8 and height 0.3 m. The 
model can implement the condition of h x = h 2 . 

In calculating >(/ /; ., an initial radiation characteristic value 
is set for each zone: attenuation (absorption) coefficient for 
the volumetric zones and the degree of blackness (reflecting 
capacity) for the surface zones. The degree of blackness of 
the brickwork surface zones (roof and walls) was set equal to 
0.8. The degree of blackness of the glass melt surface was set 
in the range of 0.8 - 0.4, decreasing toward the working 
zone. The radiation absorption coefficient was calculated us¬ 
ing the Gurvich - Mitor fonnula: 
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where k Q is the attenuation coefficient, m~ *; p is the partial 
pressure of the gaseous medium component, Pa; S ef is the ef¬ 
fective ray length, m; T is the temperature of the gas (zone), 
K; V and F are the total volume and surface area of the ga¬ 
seous zones, m 3 and m 2 , respectively. 

At the second stage the calculation of the optical-geomet¬ 
ric parameters of the zone system elements is completed. The 


matrix of the coefficients \|/ -, the radiation characteristics of 
the zones, and data on their mutual position are used to deter¬ 
mine the matrices of resolvent generalized angular radiation 
coefficients by means of solving a system of linear equa¬ 
tions of the zonal resolvent method, and then the matrix of 
radiation exchange coefficients a /: is calculated [1]. 

At the third stage a database is compiled reflecting the 
real operating conditions of the furnace, which makes it pos¬ 
sible to form a system of nonlinear equations of the resolvent 
zonal method [1], Using the restricted jets regularities, the 
aerodynamic contours and velocities of the flame near the 
enclosing brickwork and melt surfaces are determined. The 
aerodynamic contours of direct gas flows (the sequence of 
passing the zones) and the degree of combustion of fuel are 
the source data for calculating the mass exchange between 
the volumetric zones, the combustion process, and the con¬ 
vection exchange coefficients. The current value of the rela¬ 
tive air inflow to the combustion zone a and the degree of 
fuel combustion y are set in the fonn of the following func¬ 
tions [2]: 


a = a[ 1 - exp (- mnx )]; 

X = 1 - exp (-rax 2 ), 

where a is the coefficient of the air flow rate supplied for 
combustion; a = 1.1 (BT-1 and ZT-1), a = 1 (KT-1); m, n are 
the coefficients depending on the flame length; x is the dis¬ 
tance from the burner exit section, m. 

The particular values of the coefficients m and n are 
specified based on the following reasoning. It is believed that 
85% of the fuel is burnt within the limits of the inflow length 
l- m limited by the value a = 1. The rest of the fuel is con¬ 
sumed in the afterburning zone. The total length of the flame 
is l f = l-43/ jn and the integral chemical underburning of the 
fuel is equal to q 2 = 1 -% = 0.02 [2], An adaptation of the 
model is implemented for l m = 9.534 m, which corresponds 
to the extent of the fuel combustion zone approximately 
equal to the length of the working space of the furnace: 
1 { = 13.634 m. 

For all enclosing brickwork elements the convective heat 
transfer coefficient a c = 8 W/(m 2 • K). For the glass melt 
surface zones this coefficient is taken to be constant across 
the width of the calculation zone and variable along the fur¬ 
nace length: a c = 15, 15, 12, 10, and 8 W/(m 2 • K). The cal¬ 
culations included the following values of the heat capacity 
of combustion products, air, and fuel: 1.613, 1.425, and 
1.571 kJ/(m 3 • K), respectively. The fuel combustion compu¬ 
tations determined not only the composition and quantity of 
the combustion products and the radiation attenuation coeffi¬ 
cients, but also the energy parameters of the gas flow. A cor¬ 
rect determination of the volumetric density of heat emis¬ 
sions in the gaseous zones depends on the accuracy of the 
calculation of heat losses via the enclosing brickwork and 
heat assimilated in the tank. 




























Mathematical Model of a Glass-Melting Furnace with Horseshoe-Shaped Flame 


325 


In the general case the thermal flow via the enclosing 
surface j is found from the formulas: 

<lj = K j( T brj- T am); 

( « c I s ! 1 

K= Xf + — , 

^/=i A,,- a br y ^ 

where q- is the heat flow via the surface zone j, W/m 2 ; K. is 
the heat transfer coefficient via the surface zone j, 
W/(m 2 • K); T blJ and T am are the temperatures of the inner 
brickwork surface of the zone j and the ambient medium, K; 
n is the number of brickwork layers; S i is the thickness of the 
rth brickwork layer, m; X j is the thermal conductivity of the 
ith brickwork layer, W/(m 2 • K); a br; is the heat transfer co¬ 
efficient on the outer brickwork surface, W/(m 2 • K). 

For surface zones of type III the heat losses via the en¬ 
closing brickwork are taken into account by specifying the 
temperature of the ambient medium f = 30°C and the heat 
transfer coefficient. 

Let us consider in more detail the formation of the data¬ 
base for the calculation of K.. All structural elements of the 
furnace (the analog) are designed using multilayer refractory 
insulation. This ensures sufficiently low thermal losses into 
the ambient medium. The design of a thermally insulated 
dinas-brick roof and the calculation of heat losses are speci¬ 
fied in [3], It has been demonstrated that for the inner brick¬ 
work surface temperature t hr = 1580°C and t am = 30°C the ef¬ 
fectiveness of the thermal insulation of the roof 400 mm 
thick (refractory grade DSO) is characterized by the specific 
heat flows q ex = 1437.7 W/m 2 and the temperature of the ex¬ 
terior brickwork surface t ex = 112.2°C. For a noninsulated 
roof q ex = 5229.7 W/m 2 and f ex = 235.7°C. It is known that 
the roof sites adjacent to the temperatures joints and 
thermocouples are not insulated for safety reasons. This cir¬ 
cumstance is taken into account in determining the average 
value of the heat transfer coefficient for the surface zones of 
the roof taken as K = 1.29986 W/(m 2 • K). 

In determining K. the heat losses of the lateral walls of 
the flame space were taken into account not only as thermal 
conductance, but also as radiation through the inlet arches, 
loading hoppers, and peepholes. For this elements of the 
brickwork K = 4.65036 W/(m 2 • K). It significantly exceeds 
the heat transfer coefficient that would have been accepted, if 
only thermal losses via conductance had been taken into ac¬ 
count. For instance, the heat flow through insulated lateral 
walls made of badelleyite-corundum plates (BK-33) of thick¬ 
ness 200 mm is equal to q cx = 1219.5 W/m 2 . The estimated 
temperature of the external surface of the brickwork insula¬ 
tion is ? ex = 112.6°C. For comparison, note that for non-insu- 
lated walls q ex = 23,837.7 W/m 2 and t ex = 527.3°C [4], 

Special attention has been paid to the calculation of the 
heat transfer coefficient for the walls and bottom of the melt¬ 
ing tank. The lateral walls of the tank are made of palisade 


bars of thickness 250 mm (BK-33, BK-37, and BK-41). The 
upper part of the bar (200 mm) is not insulated and is sub¬ 
jected to forced air cooling. The remainder part of the bar is 
arbitrarily split in two parts. The upper part (counting from 
the melt surface) has intermittent insulation (the bar joints 
are not insulated) and the lower part has continuous insula¬ 
tion. The cooling intensity of the non-insulated part calcu¬ 
lated based on the condition t bi = 1200 - 1500°C and 
t ex = 100°C is characterized by the heat transfer coefficient 
on the exterior surface of the bar equal to a ex = 272.9 - 
365.5 W/(m 2 • K). In this case the specific heat loss is 
q ex = 19,114.3 -25,593.1 W/m 2 [5], For intermittent insula¬ 
tion q ex = 1204.8 W/m 2 and t ex = 111.9°C. For bare joints 
q ex = 16,541.2 W/m 2 and t ex = 449.5°C. A continuous insu¬ 
lation of the bar ensures q ex = 944.7 W/m 2 and f ex = 97.9°C 

[4] . The calculations of heat losses for intermittent and con¬ 
tinuous insulation are performed for the condition t br = 1450 
and 1350°C, respectively. 

The design of a heat-insulated bottom of the melting tank 
and the algorithm of its calculation are considered in detail in 

[5] . It was shown that the heat losses via the bottom in the 
melting and clarification zones can be specified as 
q ex = 1130.0 and 1205.1 W/m 2 , respectively. The specified 
heat flow values correspond to the exterior bottom tempera¬ 
tures of t ex = 119.0 and 123.6°C. The averaged values of 
heat transfer coefficients for the walls and bottom of the 
melting tank are taken equal to K = 3.571594 and 
1.02585 W/(m 2 • K), respectively. 

For the surface zones of type II representing the glass 
melt surface, the consumption items of the heat balance were 
set via the source (discharge) summands of zonal equations 
[1], This takes into account the heat consumed in glass for¬ 
mation reactions Q[, glass melt heating Q", and the compen¬ 
sation of heat losses via the brickwork of the melting tank. 
Lets us consider some specifics of the calculation of Q[ and 
Q'[. Unlike the heat balance of the furnace, the zonal model 
has to take into account not only the heat consumed in melt¬ 
ing batch, but melting cullet as well. Furthermore, the heat 
consumed in glass melt heating has to be taken into account. 
At the same time, the zonal model does not take into account 
the heat removed from the furnace with the production 
(working) glass melt flow. In this respect, the known expres¬ 
sion for calculating heat consumed in glass formation Q[ is 
transformed and has the following form. 

The heat consumed (kW) in glass formation can be cal¬ 
culated from the formula: 

Q'\ = G g \G’ b (q g + X<7u)+ G' c q c - 9 4pll ] , 

where G a is the furnace efficiency (glass melt output), 
kg/sec; G b and G c is the consumption of batch and cullet per 
1 kg of glass melt, kg/sec; q a is the summary thermal effect 
of glass formation reactions, kJ/kg; q l . = q x l + q x 2 + q x 3 is 
the heat consumption on moisture evaporation, heating the 
degassing products, and glass melting, kJ/kg of batch; q c is 
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Length coordinate x/L 



Length coordinate x/L 


Fig. 2. Temperature field on the surface of the melt (a) and roof ( b ). 
Numbers on the curves indicate glass melt temperature, °C; L is the 
furnace length equal to 13.62 m. 


the heat consumption in melting cullet, kJ/kg; g 4ph is the 
physical heat of the batch and cullet per the production of 
1 kg of glass melt, kJ/kg. 

The heat consumption of the glass-formation process de¬ 
termined from the above formula Q[ = 3381.51 kW were dis¬ 
tributed along the furnace length in the ratio of 50, 30, 15, 5, 
and 0%. The heat quantity per calculation zone is distributed 
across its width in proportion to the surface zone areas. The 
heat consumed in heating the glass melt up to 1500°C is 
equal to Q[ = 8233.912 kW. This component of heat con¬ 
sumption is distributed over all surface zones of the melting 
tank in proportion to their areas. 

The database obtained has made it possible to form a 
system of nonlinear equations of the resolvent zonal method 
[1]. Since the unknown quantities in this system are zone 
temperatures, the equations corresponding to zones of type I 
were excluded from this system. The results of solving the 
transformed system of equations are the temperatures of all 
zones II and III. The calculated temperatures are correlated 
with the temperatures field used for the calculation of the ra¬ 
diation parameters of the zones. They are taken as final only 
in the case when the difference between two temperature 
fields for all zones does not exceed the prescribed calculation 
error, whose relative value should be not more than 10" 4 . 

The results of calculating surface zone temperatures can 
be used to form regression equations making it possible to 
construct isotherm plots for these surfaces. Below are the re¬ 
gression equations of the surfaces of the melt, the roof, and 


the right and left longitudinal and front walls of the flame 
space, respectively: 

t (x, y)~ 1132.904 + 931.278* + 67.537y - 
675.024* 2 + 4.718xy - 130.01 y 2 ; 

t (x, y) = 1399.65 + 511.682* - 279.713y - 
327.564x 2 + 52.049xy + 82.726y 2 ; 

t(x,z)= 1313.145 + 761.654* +490.176z- 
542.979* 2 -275.787xz- 1646.567z 2 ; 

t (*, z) = 1278.097 + 422.825* + 408.297z- 
208.091* 2 - 128.850xz- 1567.367z 2 ; 

t(y,z)= 1575.210 — 14.301y- 350.668z + 

24.909y 2 - 361,626yz - 398.718z 2 . 

Figure 2 illustrates the data obtained with the isotherms 
for the surface of the melt and the roof. They show that the 
obtained temperature fields quantitatively and distribu¬ 
tion-wise correlate well with the real conditions of the fur¬ 
nace with the horseshoe flame direction. The melt surface 
temperature across its width varies insignificantly. In the end 
of the batch zone ( x/L ~ 0.3) it reaches 1350°C. Next, in the 
foam zone the temperature rise is somewhat slowed. The 
maximum temperature of 1450 - 1500°C is registered at the 
quellpunkt (x/L ~ 0.7) with a subsequence decrease (by ap¬ 
proximately 50°C) toward the front wall of the furnace. The 
roof temperature typically has a more perceptible heteroge¬ 
neity across the furnace width. This is due to one-sided flame 
development and high assimilation of heat in the batch melt¬ 
ing zone. Obviously, in transferring the flame the nonuni¬ 
formity of the roof temperature across the furnace will be 
leveled. The maximum roof temperature (approximately 
1580°C) is registered in the pure glass zone and is located 
exactly above the overflow edge. 

The accumulated data on optical-geometrical and radia¬ 
tion characteristics of zones, radiation and convection ex¬ 
change coefficients, and zone temperatures make it possible 
to pass to the next stage of the calculations. 

The fourth stage is the estimate of the characteristics of 
the thermal performance of the furnace, such as specific and 
integral heat flow densities and the general heat balance of 
the furnace. 

Thus, the results of the third and the fourth stages make it 
possible to reliably estimate the efficiency of the thermal per¬ 
formance of the furnace. For instance, the total quantity of 
heat contributed to the working space of the furnace due to 
the chemical heat of fuel BQ p and the physical heat of hot air 

is equal to 25213.00 kW. The processing of the calculation 
results yielded the following heat distribution among the 
thermal balance items: 

- heat assimilated by the melt surface 11,615.42 kW 
(46.07%); 

- loss via combustion products 11,923.56 kW (47.29%); 
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Fig. 3. Distribution of relative value of the flow function in the lon¬ 
gitudinal section of the tank. 

- loss via the enclosing brickwork 1337.34 kW (5.30%); 

- chemical underburning of fuel 337.36 kW (1.34%). 

The furnace efficiency is 68.86%, and the heat utilization 

coefficient is 46.07%. 

The next stage of mathematical modeling is analyzing 
the correspondence of the external heat exchange parameters 
to the conditions for producing high-quality glass melt. The 
initial data for the tank hydrodynamics problem include the 
temperature field on the melt surface (Fig. 2 a ) and the geo¬ 
metrical parameters of the tank. Computation results make it 
possible to determine the regularities of the glass melt mo¬ 
tion inside the tank (volumetric and mass velocity fields, 
multiplicity of circulation) and its temperature field. Figure 3 
shows the glass melt flow lines in the longitudinal section of 
the tank and Fig. 4 shows its temperature field in the same 
section. 

According to the definition of the flow function [1], flow 
lines represent the mass flow rate of the medium in an ele¬ 
mentary zone across the tank height. The flow function in 
Fig. 3. is expressed in fractions of the mass flow rate of mol¬ 
ten glass in the furnace neck. Since a steady furnace regime 
is considered, the flow lines coincide with the trajectories of 
the medium macroparticles. In other words, Fig. 3 shows the 
glass melt circulation contours. It can be seen that the glass 
melt in the tank has a high circulation intensity. This is re¬ 
sponsible for an insignificant nonuniformity of the melt tem¬ 
perature field, especially in the tank volume beyond the over¬ 
flow edge. In the overwhelming majority of the tank volume 
the temperature difference is not more than 15°C. The aver¬ 
age glass melt temperature across the neck section is appro¬ 
ximately 1308°C, which corresponds to the conditions of fur¬ 
naces with a high specific output. 

The glass melt temperature fluctuations in the bottom 
layers of the tank (Fig. 5) appear realistic as well. In the 
batch melting area the temperature grows in an interval of 
1218 - 1290°C. In the foam-formation area due to increasing 
surface melt temperature and due to glass melt circulation it 
remains virtually constant at a level of 1293°C. A small tem- 



Fig. 4. Temperature distribution in the longitudinal section of the 
tank. Numbers on the curves indicate glass melt temperature, °C. 
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Fig. 5. Distribution of glass melt temperature in the bottom sites 
along the tank. 

perature decrease in front of the overflow edge (1287°C) is 
due to the heat loss via the tank bottom. The glass melt tem¬ 
perature on the surface of the overflow edge reaches 
1308- 1310°C, which corresponds to the average tempera¬ 
ture of the glass melt in the neck. 

Thus, calculation results demonstrate the ample potential 
of mathematical modeling for investigating particular fur¬ 
naces and for designing new ones. 
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